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Abstract 

We  consider  a  class  of  steady-state  semil inear  reaction-diffusion 
problems  with  non-differentiable  kinetics.  The  analytical  properties  of 
these  problems  have  received  considerable  attention  in  the  literature.  We 
take  a  first  step  in  analyzing  their  numerical  approximation.  We  present 
a  finite  element  method  and  establish  error  bounds  which  are  optimal  for 
some  of  the  problems.  In  addition,  we  also  discuss  a  finite  difference 
approach.  Numerical  experiments  for  one-  and  two-dimensional  problems  are 
reported. 
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1 .  Introduction 

The  problem  we  consider  is  that  of  an  irreversible  steady-state 
reaction  taking  place  in  a  bounded  domain  Q  in  R11  (n  -  1,2,3)  with 
smooth  boundary  3Q.  The  reason  a  steady-state  occurs  is  that  the  reac¬ 
tant  used  up  in  the  reaction  in  fl  is  replenished  by  diffusion  from  the 
surrounding  region.  In  [1],  it  is  shown  that  the  problem  may  be  modeled 
by  a  scalar  equation  for  the  concentration  alone: 

(1.1)  Au  -  Xf(u(x)),  x  €  £2  (u  >  0  in  8) 

u  »  1  x  €  dft 

where  X  is  a  positive  constant  that  measures  the  ratio  of  reaction  to 
diffusion  rates.  The  function  f  measures  the  ratio  of  the  reaction  rate 
at  concentration  u  to  that  at  concentration  unity.  We  assume  that  the 
function  f  satisfies  the  following  conditions: 

(1.2)  f(t)  -  0  if  t  £  0,  f (1 )  -  1 

(1.3)  f(t)  -  tpf0(t)  for  0  (  t  <  %  0  <  m  <,  fQ(t)  <,  M  <  * 

|f0"(t)|  *  K*  0  <  p  <  1 

(1.4)  f'(t)  +  >0  for  0  <  t  £  1 . 

For  the  case  of  a  non-isothermal  reaction,  we  may  obtain  the  tempera¬ 
ture  v  of  the  reactant  from  the  relation 

(1.5)  v  -  3(1 -u )  +  1 

where  8  >  0  if  the  reaction  is  exothermic  and  8  <  0  if  it  is  endo¬ 
thermic.  We  will  primarily  be  interested  in  the  isothermal  case  (8  -  0), 


when  v  remains  identically  unity  [u  and  v  are  nondimensional 
variables]. 

The  problem  (1.1)  has  been  discussed  in  an  analytical  framework  in  [1] 
and  [3].  An  interesting  feature  analyzed  in  these  papers  is  the  exis¬ 
tence  of  a  "dead  core,"  i.e.,  a  region  in  ft  where  u  identically  van¬ 
ishes,  for  p  <1  and  X  large  enough.  Physically,  this  means  that  the 
rate  of  reaction  is  so  high  that  the  reactant  is  being  consumed  in  the 
dead  core  faster  than  it  can  be  replaced  through  diffusion  across  the 
boundary  3ft. 

The  goal  of  this  paper  is  to  take  a  first  step  in  approximating  (1.1) 
numerically.  We  will  analyze  numerical  schemes  only  for  the  case  of  pth 
order  Isothermal  reactions,  for  the  case  0  <  p  <  1 ,  where  f  is 
explicitly  given  by: 

(1.6)  f(t)  -  tp  for  t  >  0 

-  0  for  t  $  0. 

It  can  be  seen  that  this  function  satisfies  (1.2)  -  (I.1!) 

For  p  £  1 ,  this  nonlinear  problem  may  be  approximated  by  several 
methods  (for  e.g.  [2],  [4]).  However,  when  p  <  1,  f  is  not  differen¬ 
tiable  at  the  origin  (and  if  p  -  0,  it  is  not  even  continuous  there). 
Hence,  the  convergence  theorems  in  the  above  papers  fail  and  a  different 
analysis  is  required  for  the  case  of  non-differentiable  kinetics. 

The  plan  of  the  paper  is  as  follows.  In  Section  2,  we  list  existence, 
uniqueness  and  regularity  results  from  [1]  and  [3]  that  will  be  needed. 

In  Section  3,  we  present  a  finite  element  method  for  (1.1)  and  obtain  an 
error  estimate  for  its  convergence.  In  Section  4,  we  discuss  a  finite 


difference  scheme  for  the  problem.  Section  5  contains  the  results  of 
numerical  experiments  using  our  schemes. 

2.  Existence,  Uniqueness  and  Regularity  Theorems 

The  following  theorem  is  proved  in  [3]. 

2  |  (.'  . 

Theorem  2.1.  Let  f  satisfy  (1.2)  -  (1.4).  Let  ft  have  a  C 
boundary  30.  Then  there  exists  a  unique  solution  of  (1.1)  belonging  to 
C2’a(0)  where  a  -  min(p,ag). 

Since  the  function  f  given  by  (1.6)  satisfies  (1.2)  -  (1.4),  the 
above  theorem  holds  for  this  function.  In  [1],  an  alternate  proof  of 
existence  and  uniqueness  is  provided,  where  conditions  (1.3)  -  (1.4)  are 
replaced  by  the  requirement  that  f  is  positive,  continuously 
differentiable  and  monotone  increasing  on  (0,«).  This  provides  results 
for  the  function  f  in  (1.6)  for  all  p  £  0. 

By  the  maximum  principle,  it  follows  readily  that  u  $  1  in  n. 
Moreover,  the  condition  (1.2)  ensures  that  u  >  0  in  n,  which  is  es¬ 
sential  in  terms  of  physical  considerations. 

For  a  convex  two-dimensional  domain,  the  dead-core,  if  it  exists,  will 
be  convex.  The  existence  of  the  dead  core  depends  on  the  parameter  A, 
as  expressed  in  the  following  theorem,  a  slightly  modified  version  of 
which  is  proved  in  [1]. 

Theorem  2.2.  Let  f  satisfy  (1.3).  Then  there  exists  a  number  AQ  > 
0  such  that  (1.1)  has  a  dead  core  for  all  A  >  Aq. 

If  A  is  small  enough  or  if  p  >  1 ,  then  no  dead  core  will  exist. 

For  the  case  of  a  p£A  order  Isothermal  reaction  in  a  slab,  where  £2  - 
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C-d,d]  (d  is  the  width  of  the  slab)  and  f  is  given  by  (1.6),  the  dead 
core  appears  only  if 

(2.1)  vT  >  where  A  -  ^ ^ (0  £  p  <  1). 

u  1  *“P 

Moreover,  the  size  of  the  dead  core  D0  -  [-Y.Y]  is  determined  by 

(2.2)  y  -  d  -  —  . 

/T 

Similarly,  for  a  pHL  order  isothermal  reaction  in  a  ball,  there  is  a  crit¬ 
ical  value  Aq  such  that  a  dead  core  appears  if  and  only  if  A  £  AQ.  The 
above  facts  are  proved  in  [1],  where  several  other  theorems  for  the  exis¬ 
tence  and  non-existence  of  a  dead  core  are  provided. 

Finally,  it  is  shown  in  [3]  that  asymptotically,  as  A  -*•  »,  the  boun¬ 
dary  of  the  dead  core  approaches  a  smooth  surface  parallel  to  3Q  at  dis¬ 
tance  —  +  0(y) ,  6  constant. 

/T 

u  in  (1.1),  (1.6). 


p  £  (0,1) 


3.  A  Finite  Element  Approximation 


We  know  that  (3-D  has  a  unique  solution  v  €  H2(fi),  0  £  v  £  1  where 

the  region  G  -  {x  |  v(x)  -  1}  corresponds  to  the  dead  core.  The 
equivalent  weak  formulation  of  (3-1)  is  given  by: 

Find  v  €  Hg(fl)  satisfying 

(3-3)  (Vv.Vw)  -  A(g(v),w)  for  all  w  €  Hq(£1) 

p 

where  (*,♦)  denotes  the  usual  L  (fl)  inner  product. 

u  4 

Let  now  Sn  be  a  finite  dimensional  subspace  of  Hq(0)  with  the 
property 

( 3 . ^ )  inf  C |v-x|0  +  h I V-X 1 1 )  <  Chs|v|g  for  1  i  s  £  r  +  1. 

x«sh 

Here  r  denotes  the  degree  of  piecewise  polynomials  used.  Since  our 

p 

solution  is  only  known  to  lie  in  H  (Q),  we  restrict  attention  to  the 
case  r  -  1  so  that  Sh  consists  simply  of  piecewise  linear  functions. 

We  say  that  vh  6  Sh  is  an  approximate  solution  of  (3*1)  if 

(3-5)  (Wh,Vwh)  -  A(g(vh),wh)  Vwh  €  Sh. 

Theorem  3.1-  Let  g  be  given  by  (3.2).  Then  there  exists  a  func¬ 
tion  vh  €  Sh  satisfying  (3.5). 

The  proof  of  Theorem  3-1  requires  the  following  lemmas: 

Lemma  3.1.  Let  be  a  finite  dimensional  space  and  let  T  :  Mn  ■+ 

Mn  be  continuous.  Suppose  there  exists  a  sphere  Sp  with  radius  p  and 
center  at  the  origin  such  that  (Tx,x)  >0  for  x  on  Sp. 
exists  an  xQ  such  that  TxQ  «  0  and  |Xq|  £  p  . 


Then  there 


Proof ;  See,  for  example,  [7]. 

Lemma  3.2.  Let  g  be  given  by  (3.2).  Then  for  any  v1,v2  €  R 

(3.6)  |g(v-,  >  -  g(v2)|  i  |v1  -  v2|p. 

Proof.  We  may  assume  without  loss  of  generality  that  v1  >  v2. 
V1  "  v2*  (3.6)  holds  trivially.]  If  v1  >  1 ,  then  we  see  that 

U^)  -  g(v2)j  £  |g(v2)  |  i  |1  -  v2|P  $  |v1  -  v2|p. 

If  v1  <  1 ,  then 

0  i  1  -  v1  <  1  -  v2 

1-y1 

so  that  0  £  -  *  a  <  1 . 

2 

Hence,  1  -  ap  i  1  -  a  $  (1-a)p,  since  0  <  p  <  1.  This  gives 

(1-v2)p  -  ( 1 — v ^ ) p  i  (v1-v2)p 

and  (3.6)  is  proven. 

Lemma  3.3.  Let  0  <  p  <  1.  Then  for  v  >  0,  v,w  €  L2 ( £1 ) , 

(3.7)  |(vp,w)|  S  C |v |q |w |q . 

Proof.  Using  the  Schwarz  inequality,  we  have 

|(vp,w)|  i  c(J  (vP)2dx)/j  |w|0. 

Now  by  Holder’s  inequality  with  q  =  —  and  —  +  —  =  1, 
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{  v^dx  <  (/(v2p)qdx)  (J  1rdx) 

Q  “  fl  n 

<  c(J  v2dx)1/q 
n 

from  which  (3*7)  follows. 

Proof  of  Theorem  3«1« 

Define  an  operator  T  :  Sh  -*■  Sh  as  follows : 

(Tvh,wh)  -  (Wh,Vwh)  -  A(g(vh),wh)  for  all  wh  €  Sh. 

Clearly  T  is  continuous.  Moreover,  using  (3-7), 

(Tvh,vh)  >  |Wh|2  -  CA 1 1  - 

>  C(|vj|2  -  »[AP.|v’j|P]|v^|  ) 

where  A  -  measure  of  Q.  Hence,  for  |v  |1  “  P  sufficiently  large,  we 
have 

(Tvh,vh)  £  0. 

Applying  Lemma  3-1  yields  the  theorem. 

Theorem  3 .2.  (3-5)  has  a  unique  solution.  Moreover,  for  Sh  consist¬ 

ing  of  piecewise  linear  functions,  the  solution  v*1  £  0  on  P. 

Proof.  We  first  prove  uniqueness,  let  v1^1  and  vjp  be  two  solutions 
of  (3-5)  and  let  wh  -  vl}  -  v^.  Then,  from  (3-5)  we  have: 
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(Vwh,Vwh)  -  X(g(v!|)-g(Vj),  v!j-v*). 

Since  g(v)  is  decreasing  in  v,  it  follows  that 

l^nlo  <  °- 

Since  wh  -  0  on  30,  this  implies  wh  *  0,  that  is  =  v!}. 

Let  now  vh,  the  solution  of  (3.5)  satisfy 

vh(Ni)  <  0 

for  some  nodal  point  Nj  in  the  grid  on  0.  Then  there  exists  a  nodal 
point  Nq  such  that 

(3.8a)  vh(NQ)  i  vh(x)  for  all  x  £  0. 

(3.8b)  vh(N0)  <  0 

(3.8c)  v*1(Nq)  <  vh(Nj)  for  at  least  one  node  Nj  adjacent  to  Nq 

(since  vh  =  0  on  30) . 

Now,  let  <|)q  be  the  linear  basis  function  that  is  1  at  the  point  Nq 
and  0  at  all  other  nodal  points,  and  let  0q  be  its  support.  Then,  by 
(3.8a),  Vvh  is  of  the  opposite  sign  as  on  0q  and  by  (3.8c),  7vh 

is  not  Identically  zero  on  0q.  Hence,  substituting  w11  =  t|>Q  in  (3-5), 
we  obtain 

(Wh,7^)  =  j  7vh  •  V\J/Qdx  =  A(g(vh),i{;0). 

Q0 

The  left  side  is  strictly  negative  while  the  right  side  is  non-negative,  a 
contradiction. 


<7 Urt* 
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Remark.  The  non-negativity  proof  above  is  essentially  the  discrete 
maximum  principle  (using  linear  functions)  for  the  non-linear  problem 
(3-1)  with  non-negative  forcing  function  g. 

y. 

We  now  deal  with  the  question  of  convergence  of  vn  to  v. 

Theorem  3.3.  Let  v^  be  the  solution  of  (3.5)  and  v  the  solution 
of  (3-3).  Then 

|v  -  v11^  <  Ch  if  j  <  p  <  1 

<  Ch2p  if  0  <  p  <  j 


8 


s 


where  the  constant  C  is  independent  of  h  but  depends  on  u,  A,  p 
and  tl. 

Proof.  Let  x  be  as  in  (3* 4)  and  wh  (  Sh.  Then,  by  (3-3).  (3.5), 
we  have 

(3-9)  ( V(vh~x) , Vwh)  +  (g( x)-g(vh) *wh) 

*  (V(v-x) , Vwh)  +  (g( x)-g(v) ,wh) . 

Taking  wh  =  vh  -  x  in  (3*9)  gives 


(3-10) 


h  »  1 2 


|V(v“-x)|q  +  (g( x)“g(vh) ,vh-x) 


( V(v-x) . V(vh~x) )  +  (g( x)"g(v) ,vh-x) 


Since  g(v)  is  decreasing  for  v  <  1  and  zero  for  v  >  1 ,  it  follows 


(g( x)“g(vh) »vh~x)  >  0. 


Hence,  we  may  obtain  from  (3*10): 


||V( yh-x)  |q  <  >V(v-x)l0|V(vh-x)|0  +  (g  (x)-g(v)  ,vh-x) . 


Using  (3.6),  (3.7)  and  the  Poincare  inequality  gives 


|vh-X|1  <  |v— X|  1  +  CA  |v  — X|q)  • 


Using  the  approximation  property  (3.*0  together  with  the  triangle 


inequality  yields 


(3.1D 


|v-vh|,  <  C(h|v|2  ♦  Ah2p|v|P) 


from  which  the  assertion  of  the  theorem  follows. 


Remark .  The  estimate  given  in  Theorem  3-3  for  the  case  p  <  '£  is 
pessimistic.  The  numerical  experiments  in  Section  5  suggest  that  one 
gets  0(h)  convergence  for  any  p. 


A  Finite  Difference  Approach 

In  this  section,  we  present  a  finite  difference  scheme  for  (1.1), 

(1.6)  and  analyze  its  numerical  properties. 

For  the  sake  of  our  discussion,  we  take  SI  to  be  the  unit  square  0  £ 
x  (  1,  (M  y  *  1. 

Let  nh  be  a  uniform  n  «  n  finite  difference  mesh  with  mesh  spac¬ 
ing  h  and  boundary  3Ph.  Uj  j  will  indicate  an  approximate  value  of 
the  solution  of  (1.1),  (1.6)  at  (x  t ,  y  j )  =  (ih,jh)  for  0  H,j  (  N. 

uh  will  denote  the  vector  with  components  Uifj,  1  <  ^  N  1  ’ 

listed  row  by  row. 


.n'^/uVUv  v; 
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Let  be  the  usual  five-point  discrete  Laplacian  and  consider  the 

scheme 

( 4.1 )  Ahah  -  xag  on  flh,  ah  -  1  on  3tih. 

Let  m  -  (N-1)2  and  let  if1*  -  {X  |  0  £  Xj,  i  -  1 . m}.  Then  (4.1) 

may  be  expressed  as  a  non-linear  system  of  m  equations  in  m  unknowns: 

(4.2)  F(ah)  -  {f±j (ah) }  -  b 

where  F^  is  defined  on  flF1*  for  2  $  l,j  £  N-2  by 

(^.3)  Fij(ah)  *  “  ui+i,j  "  ui-i,j  “  ui,J+1  "  ui,j-1  *  4ui ,  j  *  Xh2u?,j- 

For  other  values  of  i,j,  Fjj  is  defined  once  again  by  (4.3)  except  that 
the  contributions  from  u^j  on  38h  are  taken  on  the  right  hand  side  to 
comprise  the  vector  b.  Notice  that  all  entries  in  b  will  be  non¬ 
negative.  In  what  follows,  we  will  assume  for  simplicity  that  the  compo¬ 
nents  of  F  and  0h  are  given  by  { F^ }  and  {u^  respectively,  i  - 
1,2,...,m.  We  will  be  using  the  terminology  from  [5]  and  [6]  to  define 
terms  like  inverse  isotone,  off-diagonally  antitone  and  M-function. 

The  continuous  mapping  F  :  Rm+  -»  R111  turns  out  to  be  an  M-function, 
yielding  the  following  result. 

Theorem  4.1.  The  difference  scheme  (4.1)  has  a  unique  solution  0  £ 

£  1 .  Moreover,  let  Q^.Q^  (  fiF+  be  vectors  with  all  components  0 
and  1,  respectively.  For  any  u>  €  (0,1]  define  the  (SOR)  iterates 
dk  by 


k+1  k+1  k  k 

Fi(u1  . ui-rui,ui+1 . um^  “  bi  for  ui* 


(4.4) 


ui  +  ^ : 


i  -  1 ,2, . . . ,m,  k  -  0, 1 , 


Then  the  iterates  {Qk)  starting  from  D®  and  Ojj  are  uniquely  defined, 


non-negative  and  converge  monotonically  to  Qh  from  below  and  above, 


respectively . 


Proof.  We  first  show  that  F  is  an  M-f unction.  By  Def.  2.8  of  [£], 


we  must  show  that  F  is  inverse  isotone  and  off-diagonally  antitone. 


Define  4  :  ■»  rf"  by 


(4.5) 


Mah)  *  Xh^Uj, 


1  i  i  i  m. 


Then  it  is  immediately  seen  that  $  :  ♦  hF1  is  a  continuous,  isotone, 


diagonal  mapping.  Moreover,  (4.3),  (4.5)  may  be  used  to  obtain  the 


following  splitting  for  F: 


(4.6) 


F(ah)  -  Aah  +  *(ah). 


The  m  x  m  matrix  A  is  irreducibly  diagonally  dominant  with  a^  j  $ 
0  for  i  f  j  and  a^j  >  0  for  i  =  1, — ,n.  From  this,  it  follows 


that  F  if  off-diagonally  antitone.  In  addition,  by  2.4.14,  p.  55  of 


[5],  A  is  an  M-matrix.  The  fact  that  F  is  inverse  isotone  now  follows 


easily  by  the  proof  of  13-5.6,  p.  467  of  [5]. 


By  Theorem  2.10  of  [6]  ,  F  is  strictly  diagonally  isotone.  We  now 


note  that  taking  Xq,  7q  £  *^n+  consisting  of  all  0's  and  all  1 's, 


respectively,  yields 


*0  S  ?0 


f(x0)  £  b  s  F(y0) 

and 

J  -  {x  €  R0  |  x°  £  x  £  y0}  c  tf"+. 

I 

! 

By  Theorem  3-1  of  [6],  we  find  that  for  any  <o  (  (0,1],  the  Gauss-Seidej. 

iterates  (yk)  and  {Xk)  given  by  (4.4)  and  starting  from  y0  and 

*  # 

Xq,  respectively,  are  uniquely  defined  and  satisfy  Xk  t  X  ,  yk  +  y  , 

K  ft  ~ 

F(X  )  -  F(y  )  -  b,  where  the  monotonicity  of  the  convergence  insures  0  £ 

x*,y*  i  1,  i  -  1,2 . m.  Finally,  the  fact  that  X*  -  y*  -  Qh  follows 

from  the  inverse  isotonicity  of  F,  so  that  Qh  is  unique  and  0  £  dh  £  1 . 

Remark :  It  is  obvious  that  Theorem  4.1  will  be  true  for  any  sub  and 
super  solutions  Q®  and  a®,  respectively,  for  (4.2). 

(4.4)  therefore  provides  us  with  an  Iterative  scheme  to  solve  (4.1) 
which  has  been  used  successfully  by  us  in  computations.  It  is  of  interest 
to  note  that  the  solution  dh  of  (4.1)  does  not  have  a  discrete  dead  core 
in  the  sense  that  Uj  j  i<  0  at  any  point  of  the  mesh.  For,  if  u^  j  «  0 
at  a  mesh  point,  then  by  (4.3)  we  obtain 

ui-1.j  +  ui*1.J  +  ui,j-1  +  ui,j*1  "  0 
so  that  by  non-negativity  of  the  solution, 

ui -1 , J  *  ui*1,J  "  UI , j -1  *  ui ,  j  +1  '  °* 

This  in  turn  implies  that  Dh  -  0  at  all  interior  mesh  points,  a  conse¬ 
quence  that  obviously  contradicts  (4.1)  adjacent  to 

Numerical  experiments  show  that  the  above  difference  scheme  yields 
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0(h)  convergence  in  the  discrete  H1  norm  of  the  error.  As  for  many 


difference  schemes,  a  higher  convergence  rate  is  obtained  if  it  is 


measured  at  the  nodes  alone.  The  following  difference  scheme  for  the  one¬ 


dimensional  problem  is  interesting  in  this  respect. 


(4.5) 


ui+1  -  2ut  ♦  Ui_, 


2  u?+1  +  lOu1?  +  . 

ih2  (-Ail - ^ - Li). 


For  the  case  when  A  -  12  and  p  -  %,  (4.5)  reproduces  the  true 


solution  u  -  x  exactly  at  the  mesh  points,  so  that  the  error  is  zero  at 


these  points. 


5.  Numerical  Experiments 


In  this  section  we  look  at  results  obtained  using  the  finite  element 


method.  In  the  one-dimensional  case,  the  width  of  the  dead  core  for 


(1.1),  (1.6)  can  be  analytically  determined  by  (2.1),  (2.2).  Moreover, 


when  a  -  [-1 ,+1 ],  p  -  0.5  and  A  «  12,  the  exact  solution  is  given  by 


u  -  x  ,  which  has  a  one-point  dead  core  at  the  origin.  In  this  case,  the 


error  in  the  computations  below  can  be  measured  exactly.  In  other  cases, 


the  true  solution  is  replaced  by  an  approximation  using  a  sufficiently 


small  mesh  size  [number  of  sub-intervals,  N  »  64].  The  mesh  size  h  is 


given  by  2/N. 
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p 

X 

N 

H*  error 

0.5 

12 

4 

1 .01138 

8 

0.5H50 

12 

0.3610 

16 

0.2735 

0.5 

20 

1 .  H1404 

8 

0.7923 

12 

0.5387 

16 

0.3957 

0.1 

20 

H 

1 .91129 

8 

1 .0071 

12 

0.6827 

16 

0.5035 

The  above  experiments  show  the  errors  obtained  with  linear  elements. 
It  is  observed  that  the  rate  of  convergence  in  the  H1  norm  is  0(h)  in 
all  cases.  This  corresponds  with  the  estimates  for  p  £  0.5  in  Theorem 
3-3  but  exceeds  the  rate  of  0(h2p)  predicted  for  p  <  0.5. 

The  two-dimensional  results  show  exactly  the  same  orders  of  conver¬ 
gence.  In  this  case,  experiments  are  first  performed  over  a  square  using 
linear  piecewise  polynomials  on  a  uniform  triangular  mesh.  The  next  set 
of  experiments  Involve  bilinear  functions  on  a  uniform  rectangular  mesh. 
Both  these  finite  element  spaces  satisfy  (3.1<).  N  now  represents  the 
number  of  sub-divisions  on  each  side  of  the  square  and  the  mesh  size  h 
is  once  again  2/N.  Some  sample  results  are  presented  below. 


Bilinear  Functions 


p 

A 

N 

H1  error 

0.5 

4 

4 

0.4022 

8 

0.2030 

12 

0.1352 

16 

0.0995 

24 

4 

1.7751 

8 

0.9312 

12 

0.6149 

16 

0.4596 

48 

4 

2.8526 

8 

1 .6228 

12 

1 .0570 

16 

0.8051 

0.1 

4 

4 

0.4382 

8 

0.2200 

12 

0.1468 

16 

0.1077 

10 

4 

1.1021 

8 

0.5718 

12 

0.4046 

16 

0.3145 

20 

4 

1 .9706 

8 

0.9848 

12 

0.6559 

16 

0.4815 

Remarks ; 

(a)  Dependence  on  A .  It  is  observed  that  for  large  X,  the 
increase  in  the  error  for  the  same  p  and  N  is  proportional  to  /X,  a 
predicted  by  (3-11)- 

2  ? 

(b)  L  errors ■  For  some  of  our  calculations,  the  \r  errors  con- 

2 

verged  at  the  rate  of  0(h).  In  other  cases,  the  convergence  rate  was 
lower.  We  have  not  yet  determined  conclusively  either  by  means  of  compu 
tatlon  or  analysis  what  the  correct  rate  should  be. 


(c)  The  case  p  =  0.  Although  Theorem  3-3  does  not  predict  conver¬ 
gence  in  this  case,  computationally,  we  once  again  observed  0(h)  conver¬ 
gence  in  experiments  similar  to  the  above. 

(d)  Finite  Difference  Calculations.  Computations  based  on  the  dif¬ 
ference  scheme  (4.1)  are  not  reproduced  here.  They  showed  similar  magni¬ 
tudes  of  error  and  identical  rates  of  convergence  as  the  finite  element 
case,  both  in  the  one-  and  two-dimensional  case. 

(e)  Boundedness  of  Approximate  Solutions.  In  all  our  calculations, 
we  observed  that  the  approximate  solution  always  satisfied  0  £  u^  i  1 . 

For  the  finite  difference  case,  this  is  justified  by  Theorem  4.1.  For  the 
finite  element  case,  the  discrete  maximum  principle  justifies  uh  £  1 ,  as 
shown  in  Theorem  3.2.  We  have  not,  however,  been  able  to  prove  that 

uh  *  °* 

(f)  Existence  of  the  Dead  Core.  In  the  one-dimensional  case,  it  was 
observed  computationally  that  Theorem  2.2  was  satisfied  with  AQ  =  1  2 
when  p  =  0.5.  Similar  critical  values  of  A  were  observed  for  the  two- 
dimensional  case.  In  the  two-dimensional  examples  presented,  for  A  =  4 
there  was  no  dead  core,  while  for  other  values  of  A,  a  dead  core  was 


observed . 
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